/*********************************************************************************** 
 *
 *                     ** VELSIM **
 *
 * Gaussian simulation of time-depth conversion in layer-cake model
 *
 * Inspired by work bij G.Diephuis et al. for NAM B.V.(1995-2000)
 *
 * v1.0		21/8/2011 M.A.van Ravesteijn(Raafcom B.V.)
 *			Initial version
 *
 * Format Input file (stdin)
 * -- Project name
 * -- Grids file name (Grid format x,y,twt1,twt2....twtnrlayers, other grid parameters 
 *    defined in following rows)
 * -- Gridvvn Grids file array index (layer 1..10) (value -1 no grid parameter)
 * -- Gridvvndev Grids file array index (layer 1..10) (value -1 no grid parameter)
 * -- Gridtwtdev Grids file array index (layer 1..10) (value -1 no grid parameter)
 * -- Avz(Average depth) results filename
 * -- Sdz(Standard Deviation depth) results filename
 * -- Report filename
 * -- nrlayers
 * -- Simulation option
 *    (1 = Single pont single run
 *     2 = Single point simulation
 *     3 = All points single run.
 *     4 = all points full simulation)
 * -- Nr of realizations in simulations(2-2000)
 * -- Nr of layers velocity model.
 * -- (Per layer)
 *   -- Layernr.
 *   -- Velocity type(V0/Vint) 1 = Vint, 2 = Reg Vint, 3 = V0 and k, 4 = Reg V0 and k
 *   -- Velocity(V0/Vint)(m/s).
 *   -- Standard Deviation Velocity(m/s).
 *   -- K(1/s)
 *   -- Standard deviation TWT(ms)
 ***********************************************************************************
 */

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <time.h>

#define CONS_V			1
#define REG_V			2
#define CONS_V_K		3
#define REG_V_CONS_K	4

FILE *fd_rep_vs,*fd_grids,*fd_avz_vs,*fd_sdz_vs;
FILE *fp_log_vs;

double random_vs();
double invnor_vs(double);


int main()
{
   double *z1_vs,*z2_vs,*vn_vs,*v1_vs,*v2_vs,*t1_vs,*t2_vs;

   float  velo_vs[10],sdvn_vs[10],rka_vs[10];
   float  sdtwt_vs[10],twt_vs[10],temp2_vs[40];
   float  temp3_vs[22],rav_vs[12],rsd_vs[12],regv_vs[10];
   float  xcoord_vs,ycoord_vs,xcoordv_vs,ycoordv_vs;
   float  fvalue_vs;
	
   double sum_vs,sum2_vs,deltaz_vs,deltat_vs;
   double diff_vs,dumt_vs,xt_vs,qt_vs,dumv_vs,xv_vs,qv_vs,tsq_vs,ssq_vs;

   int    ivm_vs[20];
   int    ln_vs[20];
   int    isimul_vs,nrsimul_vs,np_vs,nlay_vs,nsim_vs;
   int	  nrgrid;
   int    igridv0[10],igridvvndev[10],igridtwtdev[10];
   int    i,j,k,l,i1,v;
   int    iloop_vs; 
   int    nv_vs;
   int    nfirst_vs;

   char   projnam[101],fgrids[101],favz_vs[101],fsdz_vs[101],frep_vs[101];
   char   buf[256];
   char   string[500],string2[20];
   char   string3[20];
   char   *p,*q,*end;
   char   *starttime,*endtime;

   typedef int Boolean;
   enum { 	False = 0,
		    True  = 1};
   Boolean disc_vs,gridvvn[10],gridvvndev[10],gridtwtdev[10];

   disc_vs = False;   
   nv_vs = 10;        
   for (v=0;v<nv_vs;v++) {
        gridvvn[v]=False;
        gridvvndev[v]=False;
		gridtwtdev[v]=False;
   }
    srand(RAND_MAX/2);
	nlay_vs   = 0;
	isimul_vs = 0;
	xcoord_vs  = 0.0;
	ycoord_vs  = 0.0;
	xcoordv_vs = 0.0;
	ycoordv_vs = 0.0;
    fp_log_vs = fopen("velsim.log","w");
	if(	(z1_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(z2_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(vn_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(v1_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(v2_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(t1_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(t2_vs = (double *) calloc(2000,sizeof(double))) == NULL 
	   )
	{
		fprintf(fp_log_vs,"Not possible to allocate memory for arrays\n");
		exit(1);
	}
   
   fprintf(fp_log_vs,"In velsim.c\n");
   char startdateStr [9];
   char starttimeStr [9];
   _strdate( startdateStr);
   fprintf(fp_log_vs, "The current date is %s \n", startdateStr);
   _strtime( starttimeStr );
   fprintf(fp_log_vs, "The current time is %s \n", starttimeStr);
//  Read Projectname
   gets(projnam);
   projnam[strlen(projnam)]='\0';	
   fprintf(fp_log_vs,"projnam [%s]\n",projnam);
//  Read filename of Grids file
   gets(fgrids);
   fgrids[strlen(fgrids)]='\0';        
   fprintf(fp_log_vs,"fgrids[%s]\n",fgrids);
   fprintf(fp_log_vs,"nv_vs[%d]\n",nv_vs);
//  Read gridparameters
	for (v=0;v<nv_vs;v++) {
	   scanf("%d",&igridtwtdev[v]);
   }
   for (v=0;v<nv_vs;v++) {
	   scanf("%d",&igridv0[v]);
   }
   for (v=0;v<nv_vs;v++) {
	   scanf("%d",&igridvvndev[v]);
   }
/*
 * Read output filenames
 */
   scanf("\n");
   gets(favz_vs);
   favz_vs[strlen(favz_vs)]='\0';	
   fprintf(fp_log_vs,"favz_vs = %s\n",favz_vs);
     
   scanf("\n");
   gets(fsdz_vs);
   fsdz_vs[strlen(fsdz_vs)]='\0';	
   fprintf(fp_log_vs,"fsdz_vs = %s\n",fsdz_vs);
          
   gets(frep_vs);
   frep_vs[strlen(frep_vs)]='\0';	
   fprintf(fp_log_vs,"frep_vs = %s\n",frep_vs);     
      
   if((fd_rep_vs=fopen(frep_vs,"w")) == NULL) {
      fprintf(fp_log_vs,"**ERROR** REP file %s Open failed\n",frep_vs);
	  exit(1);
   }
//   Read nrlayers,simulationoption,nr.realisations,layer parameters
	scanf("%d",&nlay_vs);
	if (nlay_vs > 10) {
	    nlay_vs=10;
		fprintf(fp_log_vs," Nr.layers>10 set to 10\n");
	}   
	scanf("%d",&isimul_vs);
	scanf("%d",&nrsimul_vs);
	for (i=0;i<nlay_vs;i++) {
	   scanf("%3d%1d %f %f %f %f",&ln_vs[i],&ivm_vs[i],
		  &velo_vs[i],&sdvn_vs[i],&rka_vs[i],&sdtwt_vs[i]);
//   convert standard deviation twt from ms to s 
	   sdtwt_vs[i] = sdtwt_vs[i]/1000.0;
	}
 
	fprintf(fd_rep_vs,"**********************************************\n");
	fprintf(fd_rep_vs,"*** REPORT on VELSIM RUN *******************\n");
	fprintf(fd_rep_vs,"**********************************************\n");
	fprintf(fd_rep_vs,"\n");
	fprintf(fd_rep_vs,"*** Input ***\n");
	fprintf(fd_rep_vs,"Project : %s\n",projnam);
    fprintf(fd_rep_vs,"Start Simulation %s %s\n", startdateStr, starttimeStr);
    fprintf(fd_rep_vs,"Nr of Layers : %2d\n");
	fprintf(fd_rep_vs,"Simulation Type : %d (1 = Single Point Single Run\n"
		              "                      2 = Single Point Simulation\n"
					  "                      3 = All Points Single Run\n"
					  "                      4 = All Points Simulation\n");
	fprintf(fd_rep_vs,"Nunber of Realisations : %4d\n",nrsimul_vs);
	fprintf(fd_rep_vs,"\n");
	fprintf(fd_rep_vs,"Layer Parameters\n");
	fprintf(fd_rep_vs,
    " ----------------------------------- \n"
	" |LY| V | VINT | STDV |  K   | STDV |\n"
	" |NR|TYP|  V0  |  V   |      |  TWT |\n"
	"  -- --- ------ ------ ------ ------ \n");
	for (i=0;i<nlay_vs;i++) {
	   fprintf(fd_rep_vs,"\n |%2d|%3d|%6.1f|%6.1f|%6.4f|%6.4f|\n",
	    ln_vs[i],ivm_vs[i],velo_vs[i],sdvn_vs[i],rka_vs[i],
	    sdtwt_vs[i]);
	}
	fprintf(fd_rep_vs,"\n");
//   Open Grids file
	if((fd_grids=fopen(fgrids,"rb")) == NULL) {
	   fprintf(fp_log_vs,"**ERROR** Grids file %s Open failed\n",fgrids);
	   exit(1);
	}
//   Report Velocity types
	for (i=0;i<nlay_vs;i++) {
	   i1=i+1;
	   switch(ivm_vs[i])
           {
	     case CONS_V:
	        fprintf(fd_rep_vs,"\n Layer %2d has a Constant VINT",i1);
	        break;
	     case REG_V:
	        fprintf(fd_rep_vs,"\n Layer %2d has a Regionally varying VINT",i1);
	        break;
	     case CONS_V_K:
	        fprintf(fd_rep_vs,"\n Layer %2d has Constant V0 and K",i1);
	        break;
	     case REG_V_CONS_K:
	        fprintf(fd_rep_vs,"\n Layer %2d has Regionally varying V0 and Constant K",i1);
	        break;
	     default: 
			 fprintf(fp_log_vs,"**ERROR** Unknown velocity type %d for Layer %2d : %d\n",ivm_vs[i],i1);
	        exit(1);
	        break;
	   }
   	}

	fprintf(fd_rep_vs,"\n");
//   Check if Grid VVN is needed
	for (i=0;i<nlay_vs;i++) {
	   if( ivm_vs[i]== REG_V || 
	       ivm_vs[i]== REG_V_CONS_K )
	   { 
	     gridvvn[i] = True;
	   }
	}
//   Check if VVN Dev Grid file is needed
    fprintf(fd_rep_vs,"\n");
	for (i=0;i<nlay_vs;i++) {
	   i1=i+1;
	   if( igridvvndev[i] > -1 ) {
	      fprintf(fd_rep_vs,"\n Layer %2d has Regionally varying V0/VINT Uncertainty",i1);
 	      gridvvndev[i] = True;
 	   }
 	   else {
 	      fprintf(fd_rep_vs,"\n Layer %2d has Constant V0/VINT Uncertainty",i1);
 	   }
	}
//   Check if TWT Dev Grid file is needed
    fprintf(fd_rep_vs,"\n");
	for (i=0;i<nlay_vs;i++) {
	   i1=i+1;
	   if(igridtwtdev[i] > -1) {
	      fprintf(fd_rep_vs,"\n Layer %2d has Regionally varying TWT Deviation [-]",i1);
 	      gridtwtdev[i] = True;
 	   }
 	   else {
 	      fprintf(fd_rep_vs,"\n Layer %2d has Constant TWT Deviation [-]",i1);
 	   }
	}

	if((fd_avz_vs=fopen(favz_vs,"w")) == NULL) {
	   fprintf(fp_log_vs,"**ERROR** AVZ file %s Open failed\n",favz_vs);
	   exit(1);
	}

	if((fd_sdz_vs=fopen(fsdz_vs,"w")) == NULL) {
	   fprintf(fp_log_vs,"**ERROR** SDZ file %s Open failed\n",fsdz_vs);
	   exit(1);
	}
//   set number of realisations following simulation type

	if (isimul_vs == 1 || isimul_vs == 3) {
	   nsim_vs = 1;
	}
	else if (isimul_vs == 2 || isimul_vs == 4) {
	   nsim_vs = nrsimul_vs;
	}
//   set standard deviations to 0 if single run
	if (nsim_vs == 1) {
	   for (i=0;i<nlay_vs;i++) {
	      sdvn_vs[i] = 0.0;
	      sdtwt_vs[i]= 0.0;
	   }
	}

	fclose(fd_grids);

	if((fd_grids=fopen(fgrids,"rb")) == NULL) {
	   fprintf(fp_log_vs,"**ERROR** Grids file %s Open failed\n",fgrids);
	   exit(1);
	}

	
	fprintf(fd_rep_vs,"\n");
//   Begin loop over xy points
	np_vs =0;
	nfirst_vs = 0;
    while(fgets(string,500,fd_grids) !=0) {
	   disc_vs=False;
       np_vs++;
	   sscanf(string,"%f %f",&xcoord_vs,&ycoord_vs);
       p    = string;
       end  = string + strlen(string) - 1;
       i    = 0;
       while (*p != '\n')
       {  q    = end;
          fvalue_vs = (float) strtod(p, &q);
          if (strcmp(p,q)==0){
             *p='\n';
             break;
          }
          if (i > 1) {
             temp2_vs[i-2] = fvalue_vs ;                    
          }
          i++;
          p    = q;
	   }
//      Set xy output arrays
	   rav_vs[0] = xcoord_vs;
	   rav_vs[1] = ycoord_vs;
	   rsd_vs[0] = xcoord_vs;
	   rsd_vs[1] = ycoord_vs;

	   for (i=0;i<nlay_vs;i++) {
	      if (temp2_vs[i]<1.0 || temp2_vs[i]>9998.0) {
		 fprintf(fd_rep_vs,"Point X= %10.1f Y= %10.1f"
			 " Discarded because of invalid TWT's\n",xcoord_vs,ycoord_vs);
	      disc_vs = True;
	      }
	   }
//      Check integrity TWT's of Point
       for (i=1;i<nlay_vs;i++) {
	      diff_vs = temp2_vs[i] - temp2_vs[i-1];
	      if (diff_vs>-2.0 && diff_vs<0.0) temp2_vs[i]=temp2_vs[i-1];
	      if (diff_vs<-2.0) {
		 fprintf(fd_rep_vs,"Point X= %10.1f Y= %10.1f"
			 " Discarded because of Discrepancy in TWT's between Layers\n",xcoord_vs,ycoord_vs);
		 disc_vs = True;
	      }
	   }
//	   READ REGIONAL VELOCITIES FROM FILE IF CALLED FOR BY THE MODEL
	   for (i=0;i<nlay_vs;i++) {
	      temp3_vs[i] = 0.0;
           }
       for (v=0;v<nv_vs;v++) {      
	      if (gridvvn[v]) {
			    temp3_vs[v] = temp2_vs[nlay_vs + igridv0[v]];
          }
       }
	   for (i=0;i<nlay_vs;i++) {
		  if (gridvvn[i]) {
			  regv_vs[i] = temp3_vs[i];
			  if (regv_vs[i] < 1.0 || regv_vs[i] > 9998.0) {
                 fprintf(fd_rep_vs,"Point X=%10.1f Y=%10.1f Discarded because of Invalid Reg.V\n",
					 xcoord_vs,ycoord_vs);
		         disc_vs=True;
		      }
	      }
	   }
//      Read Velocity Uncertainty Grid Value from Grids file record    
	   for (i=0;i<nlay_vs;i++) {
	      temp3_vs[i] = 0.0;
       }
       for (v=0;v<nv_vs;v++) {      
	      if (gridvvndev[v]) {
     		 temp3_vs[v] = temp2_vs[nlay_vs + igridvvndev[v]];
          }
       }

	   for (i=0;i<nlay_vs;i++) {
		  if (gridvvndev[i] && nsim_vs>1) {
			 sdvn_vs[i] = temp3_vs[i];
		  }
		  else if (nsim_vs==1) {
			 sdvn_vs[i]=0;
          }
       }
//      Read TWT uncertaintie from Grids file           
	   for (i=0;i<nlay_vs;i++) {
	      temp3_vs[i] = 0.0;
       }
       for (v=0;v<nv_vs;v++) {      
	      if (gridtwtdev[v]) {
			  temp3_vs[v] = temp2_vs[nlay_vs + igridtwtdev[v]];
          }
       }
	   for (i=0;i<nlay_vs;i++) {
		  if (gridtwtdev[i] && nsim_vs>1) {
//         Convert sd twt from ms to s
			 sdtwt_vs[i] = temp3_vs[i]/1000.0;
          }
	      else if (nsim_vs==1) {
			 sdtwt_vs[i]=0.0;
	  	  }
      }
		 

	   if (!disc_vs) {
	      for (i=0;i<nlay_vs;i++) {              
//          Convert twt's from ms to s
		         twt_vs[i] = temp2_vs[i]/1000.0;
              
          }


	      for (j=0;j<nsim_vs;j++) {
		     v1_vs[j] = 0.0;
		     v2_vs[j] = 0.0;
		     t1_vs[j] = 0.0;
		     t2_vs[j] = 0.0;
		     z1_vs[j] = 0.0;
		     z2_vs[j] = 0.0;
	      }
//         Start loop over Layers
	      for (i=0;i<nlay_vs;i++) {
//            Copy points from previous layer
             for (j=0;j<nsim_vs;j++) {
	            t1_vs[j] = t2_vs[j];
	            v1_vs[j] = v2_vs[j];
	            z1_vs[j] = z2_vs[j];
	         }

/*
 *  CHECK WHETHER THE TWT OF THE BOTTOM OF THE CURRENT LAYER IS EQUAL
 *  TO THE TWT OF THE BOTTOM OF THE PREVIOUS LAYER, I.E. IN THE CASE
 *  OF ABSENCE OF THE LAYER :
 *  IN THAT CASE BOTH THE AVERAGE AND STANDARD DEVIATION ARE
 *  DIRECTLY COPIED FROM THE RESULTS OF THE SIMULATION OF THE BOTTOM
 *  OF THE LAYER ABOVE.
 *  IF THE DIFFERENCE IN TWT IS LESS THAN TWO TIMES THE GIVEN STANDARD DEVIATION
 *  THEN THE STANDARD DEVIATION OF RESULTS IS COPIED BUT THE AVARAGE IS ADAPTED
 *  BY USING THE INPUT PARAMETERS ; INTERNAL SIMULATION RESULTS ARE ADAPTED
 *  ACCORDINGLY.
 */

		 if (i==0) {
		    deltat_vs = twt_vs[0];
		 }
		 else {
		    deltat_vs = twt_vs[i] - twt_vs[i-1];
		 }
		 if (i>0 && deltat_vs <= (2*sdtwt_vs[i])) {
		    rsd_vs[i+2] = rsd_vs[i+1];
		    if (twt_vs[i] == twt_vs[i-1]) {
		       rav_vs[i+2] = rav_vs[i+1];
		    }
		    else {
		       switch(ivm_vs[i])
		       {

			   case CONS_V: 
			    deltaz_vs = deltat_vs * velo_vs[i] /2;
			    break;

			   case REG_V: 
			    deltaz_vs = deltat_vs * regv_vs[i] /2;
			    break;

			   case CONS_V_K: 
			    deltaz_vs = ((rav_vs[i+1]*rka_vs[i]+velo_vs[i])/rka_vs[i])*
				     (exp( (double) (rka_vs[i]*deltat_vs/2))-1);
    			sprintf(string3,"%7.1f",deltaz_vs);
	   			string3[strlen(string3)] = '\0';
				if (strcmp(string3,"    NaN") == 0) {
				   printf("NaN in deltaz_vs  rka_vs %d %7.4f\n",i,rka_vs[i]);
				}
		        break;
              
			   case REG_V_CONS_K: 
			    deltaz_vs = ((rav_vs[i+1]*rka_vs[i]+regv_vs[i])/rka_vs[i])*
				     (exp( (double) (rka_vs[i]*deltat_vs/2))-1);
			    break;

			   default:
		        break;
		       }
		       rav_vs[i+2] = rav_vs[i+1] + deltaz_vs;
		       for (j=0;j<nsim_vs;j++) {
			      t2_vs[j] = t1_vs[j] + deltat_vs;
			      v2_vs[j] = (2*deltaz_vs)/deltat_vs;
			      z2_vs[j] = z1_vs[j]+deltaz_vs;
		       }
		    }
		 }
		 else {

/*
 *  THE TWT OF THE BOTTOM OF THE CURRENT LAYERS DIFFERS WITH THE BOTTOM
 *  OF THE PREVIOUS LAYER (I.E.) ITS OWN TOP BY MORE THAN TWO TIMES THE
 *  GIVEN STANDARD DEVIATION ; THE SIMULATION PROCESS CAN START
 */

/*
 *  GENERATE ARRAY OF TIME VALUES IF STANDARD DEVIATION .NE. ZERO
 *  TIME VALUES ARE CHOSEN SUCH THAT THE SIMULATED TIMES OF THE
 *  BOTTOM CAN NEVER BE SMALLER THAN THE TIME OF THE TOP OF AN
 *  INTERVAL  
 */

		    if (sdtwt_vs[i] > 0.0001) {
		       if (i > 0) {
				  for (k=0;k<nsim_vs;k++) {
				     diff_vs = -999.9;
				     iloop_vs = 0;
				     while (diff_vs <= 0.0 && iloop_vs < 500) {
                        iloop_vs++;
					    qt_vs   = random_vs();
					    xt_vs   = invnor_vs(qt_vs);
					    dumt_vs = twt_vs[i] + xt_vs*sdtwt_vs[i];
					    diff_vs = dumt_vs - t1_vs[k];
			         }
			         t2_vs[k] = dumt_vs;
                     if (diff_vs <= 0) t2_vs[k] = t1_vs[k] + 0.0001;
			      }
		       }
		       else {
			      for (k=0;k<nsim_vs;k++) {
			         qt_vs   = random_vs();
			         xt_vs   = invnor_vs(qt_vs);
			         t2_vs[k] = twt_vs[i] + xt_vs*sdtwt_vs[i];
			   }
		    }
		    }
		    else {
		       for (k=0;k<nsim_vs;k++) {
			      t2_vs[k] = twt_vs[i];
		       }
		    }
			 /* Serge van Gessel test stuff
			 int nrt = 40;
			 int nrv = 50;
			 int k2 = 0;
			 int k3 = 0;
             int kt = 0;
			if (i > 0) {
				for (k=0;k<nrt;k++) {
				     diff_vs = -999.9;
				     iloop_vs = 0;
				     while (diff_vs <= 0.0 && iloop_vs < 500) {
                        iloop_vs++;
					    qt_vs   = random_vs();
					    xt_vs   = invnor_vs(qt_vs);
					    dumt_vs = twt_vs[i] + xt_vs*sdtwt_vs[i];
					    diff_vs = dumt_vs - t1_vs[k];
			         }
					 for (k2=0;k2<nrv;k2++) {
                        kt = (k*nrv) + k2;
						t2_vs[kt] = dumt_vs;
						if (diff_vs <= 0) t2_vs[kt] = t1_vs[kt] + 0.0001;
					 }
			      }
				}
		       else {
			      for (k=0;k<nrt;k++) {
			         qt_vs   = random_vs();
			         xt_vs   = invnor_vs(qt_vs);
					for (k2=0;k2<nrv;k2++) {
                        kt = (k*nrv) + k2;
						t2_vs[kt] = twt_vs[i] + xt_vs*sdtwt_vs[i];
						
					 }
				  }
			   }*/
/*
 *  Fill Velocity Draws
 */

			  double vel = 0.0;
			  if( ivm_vs[i]== REG_V || 
	              ivm_vs[i]== REG_V_CONS_K )
              {
				  vel = regv_vs[i];
			  }
			  else
			  {
				  vel = velo_vs[i];
			  }
			  if (sdvn_vs[i] > 0.01) {
			     for (k=0;k<nsim_vs;k++) {
    				dumv_vs = 0.0;
					while (dumv_vs <= 0.0) {
					   qv_vs = random_vs();
					   xv_vs = invnor_vs(qv_vs);
					   dumv_vs = vel + xv_vs*sdvn_vs[i];
					}
					v2_vs[k] = dumv_vs;
					vn_vs[k] = dumv_vs;
			     }
			  }
			  else {
			     for (k=0;k<nsim_vs;k++) {
					v2_vs[k] = vel;
					vn_vs[k] = vel;
			     }
			  }
			/* Serge van Gessel test stuff
			      for (k=0;k<nrv;k++) {
    				dumv_vs = 0.0;
					while (dumv_vs <= 0.0) {
					   qv_vs = random_vs();
					   xv_vs = invnor_vs(qv_vs);
					   dumv_vs = vel + xv_vs*sdvn_vs[i];
					}
				  
					for (k2=0;k2<nrt;k2++) {
                        kt = (k2*nrv) + k;
						v2_vs[kt] = dumv_vs;
						vn_vs[kt] = dumv_vs;
					 }
					}*/

/*
 *  PERFORM V0&K DEPTH CONVERSION AND CALCULATE PSEUDO INTERVAL
 *  VELOCITIES
 */
			  
			 if ((ivm_vs[i] == CONS_V_K ||
			      ivm_vs[i] == REG_V_CONS_K) &&
				  rka_vs[i] > 0.0000001) {
			    for (k=0;k<nsim_vs;k++) {
				  if (fabs(rka_vs[i])>0.0000001) {
					 z2_vs[k] = z1_vs[k] + ((vn_vs[k] + z1_vs[k]*rka_vs[i])/rka_vs[i]) *
						 (exp( (double) (rka_vs[i]*(t2_vs[k]-t1_vs[k])/2.0) )
						  - 1.0);
				  }
				  else {
					 z2_vs[k] = z1_vs[k] + (vn_vs[k] * (t2_vs[k]-t1_vs[k])/2.0);
				  }
				  v2_vs[k] = 2.0 * (z2_vs[k]-z1_vs[k])/(t2_vs[k]-t1_vs[k]);
			   }
			 }
		 
			  

/*
 *  PERFORM THE DEPTH CALCULATIONS
 */
		    for (k=0;k<nsim_vs;k++) {
		       z2_vs[k] = z1_vs[k] + v2_vs[k]*(t2_vs[k]-t1_vs[k])/2;
		    }

/*
 *  CALCULATE MEAN AND STANDARD DEVIATION OF RESULTS
 */
		    sum_vs  = 0.0;
		    sum2_vs = 0.0;
			double sumsqr = 0.0;
		    for (k=0;k<nsim_vs;k++) {
		       sum_vs  = sum_vs + z2_vs[k];
			   tsq_vs  = z2_vs[k];
		       sum2_vs = sum2_vs + tsq_vs*tsq_vs/(double) nsim_vs;                            
		    }
		    rav_vs[i+2] = sum_vs / (double) nsim_vs;
			for (k=0;k<nsim_vs;k++) {
				sumsqr = sumsqr + ((z2_vs[k] - (sum_vs/(double) nsim_vs))*(z2_vs[k] - (sum_vs/(double) nsim_vs)));
			}
			tsq_vs      = sumsqr / (double) nsim_vs;
//          tsq_vs      = sum_vs/ (double) nsim_vs;
            ssq_vs      = sum2_vs - tsq_vs*tsq_vs;
            if(ssq_vs<0.0) {
                rsd_vs[i+2]=0.0;
            }
            else {
//            rsd_vs[i+2] = sqrt(ssq_vs);
			  rsd_vs[i+2] = sqrt(tsq_vs);
            }

/*
 *  END OF THE CALCULATIONS ON ONE LAYER
 */
		 } /* end of if(i>0&& deltat_vs < .....   else  */
	      }

/*
 *  WRITE RESULTS TO FILE
 */
              
	      if(nfirst_vs != 0) {
                 fprintf(fd_avz_vs,"\n%10.1f%10.1f",rav_vs[0],rav_vs[1]);
              }
              else {
                 fprintf(fd_avz_vs,"%10.1f%10.1f",rav_vs[0],rav_vs[1]);
              }

	      for (k=0;k<nlay_vs;k++) {
			   sprintf(string3,"%8.1f",rav_vs[k+2]);
				string3[strlen(string3)] = '\0';
				if (strcmp(string3,"    NaN") == 0) {
					rav_vs[k+2]=9999.9;
				}
		      fprintf(fd_avz_vs,"%8.1f",rav_vs[k+2]);
	      }
              if (nfirst_vs != 0) {
	         fprintf(fd_sdz_vs,"\n%10.1f%10.1f",rsd_vs[0],rsd_vs[1]);
              }
              else {
	         fprintf(fd_sdz_vs,"%10.1f%10.1f",rsd_vs[0],rsd_vs[1]);
              }
			  nfirst_vs  = 1;

	      for (k=0;k<nlay_vs;k++) {
			   sprintf(string3,"%8.1f",rsd_vs[k+2]);
				string3[strlen(string3)] = '\0';
				if (strcmp(string3,"    NaN") == 0) {rsd_vs[k+2]=9999.9;}
		      fprintf(fd_sdz_vs,"%8.1f",rsd_vs[k+2]);
	      }

	   } /* end of if (!disc_vs) compound)  */
/*
 *  END OF CALCULATIONS BELONGING TO ONE (X,Y) POINT
 */
	}
/*
 *  CLOSE FILES
 */
	      free(z1_vs);
	      free(z2_vs);
	      free(vn_vs);
	      free(v1_vs);
	      free(v2_vs);
	      free(t1_vs);
	      free(t2_vs);

	      fclose(fd_grids);
              
              
          fprintf(fd_avz_vs,"\n");
	      fclose(fd_avz_vs);
	      fclose(fd_sdz_vs);
		  fprintf(fd_rep_vs,"\n THE NUMBER OF CDP'S READ IN IS : %4d",np_vs);
	      fprintf(fd_rep_vs,"\n NORMAL EXIT OF PROGRAM");
		  fprintf(fd_rep_vs, "\n Start Simulation %s %s", startdateStr, starttimeStr);
		  char enddateStr [9];
		  char endtimeStr [9];
		  _strdate( enddateStr);
  		  _strtime( endtimeStr );
		  fprintf(fd_rep_vs, "\n End Simulation   %s %s\n", enddateStr, endtimeStr);
		  
    
	      fclose(fd_rep_vs);
          fclose(fp_log_vs);
             
	      exit(0);
     }


/********************************************************************
 *       Function : invnor_vs
 *
 *       Function to compute a function value (y * stdev)
 *       given a probaility for a gaussian distrtibution
 *
 *       Derived from www.johndcook.com/normal_cdf_inverse.html 
 *       (based on Abrhamovic)
 *
 ********************************************************************/

double invnor_vs(double q)
{
   double tt,x,y;
   double c0 = 2.515517;
   double c1 = 0.802853;
   double c2 = 0.010328;
   double d1 = 1.432788;
   double d2 = 0.189269;
   double d3 = 0.001308;


   if (q <= 0.5) {
        tt = sqrt(-2.0*log(q));
   }
   else  {
        tt = sqrt(-2.0*log(1.0-q));
   }
   y = tt - ((c0 + ((c1 + c2*tt)*tt))/
	    (1 + (d1 + (d2 + d3*tt)*tt)*tt));

   if (q <= 0.5) {
      return(y);
   }
   else {
      return(-1.0*y);
   }
}

/********************************************************************
 *      PROGRAM : random_vs
 *
 *  SUBROUTINE TO GENERATE A RANDOM NUMBER BETWEEN 0.0 AND 1.0
 *  BOTH EXCLUSIVE.
 *  
 ********************************************************************/
 
double random_vs()
{
    double x;
    int    i;
    x = 0.0;
    while (x < 0.00000001 || x > 0.999999999) { 
/*  while (x < 0.00001 || x > 0.99999) { */
       x = (double) rand() / (double) RAND_MAX;
    }
    return x;
}
